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INTRODUCTION 


When  expensive,  highly  reliable  systems  are  produced  for  a 
customer,  the  customer  often  demands  more  assurance  than  a  point 
estimate  of  the  reliability  from  component  by  component  test  data. 
Many  methods  [IDS  1979]  have  been  proposed  for  computing  lower 
confidence  bounds  for  series  systems  but  few  methods  have  been 
extended  to  general  coherent  systems  including  pseudo-series 
[IDS  1980] ,  a  structure  which  may  be  used  to  model  certain  types  of 
missiles. 


The  likelihood  ratio  method,  introduced  by  Madansky  [M]  and 
extended  to  all  coherent  systems  by  Myhre  and  Saunders  [M,S],  forms 
the  basis  of  many  methods  that  have  proven  useful  for  series  systems 
[IDS  1979]  and  pseudo-series  [IDS  1980].  The  computations  of  this 
method  are  difficult  to  do  by  hand  for  series  systems  and  can  be  so 
complicated  for  other  systems  as  to  necessitate  a  computer  program. 
Subroutines  which  do  these  computations  for  pseudo-series  and  other 
types  of  coherent  systems  are  the  focus  of  this  report. 
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THE  LIKELIHOOD  RATIO  METHOD 


In  [M] ,  Madansky  introduced  the  likelihood  ratio  method  of 
computing  confidence  intervals  for  series  and  parallel  systems  with 
binomial  component  data;  this  method  is  based  on  a  result  of  Wilks 
that  -21n  L(r)  is  distributed  asympotically  as  a  chi-squared  with  one 
degree  of  freedom,  where  L(r)  is  the  likelihood  ratio: 

L(r)  =  sup  B(Xi,nn.;p.) 

sup  b(x  vrLp-p  j)  , 

where  B(^,n^;f^)  is  the  binomial  distribution  and  where  the  sup  in 
the  numerator  is  over  all  p  such  that  r=H  pj  •  while  the  sup  in  the 
denominator  is  over  all  possible  values  of  p  and  is  achieved  when  p  is 
the  maximum  likelihood  estimate  p^=x ^/n^  .  In  [M]  Madansky  shows  that 
the  numerator  can  be  made  a  function  of  5  ,  a  Lagrange  multiplier, 
with  5^  =  (X£-<5)/(n  £-<5)  and  r  =  np^.  Hence  it  is  straightforward  to 
solve  for  <5 1<0  and  5  2>0  by  finding  the  two  solutions  of 
In  L(p(6) )=-x2(1)/2  ,  where  x2(l)  is  the  1-y  percentile  of  the 

chi-sqaured  distribution  with  one  degree  of  freedom.  The  resulting 
confidence  interval  is  S(r)={r:  h(p^)<r<h(p(  5  ))}  ,  where  h(p)  is 
the  reliability  function  of  the  system  evaluated  at  the  vector  p  and 
Pi(5)  =  (xi-5)/(ni^  ). 

In  [M,S]  Myhre  and  Saunders  extended  this  method  • to  general 
coherent  systems  with  the  resulting  computations  being  more 
complicated  than  those  in  the  series  case.  Using  the  method  of 
Lagrange  with  multiplier  <5  to  find  the  sup  of  the  numerator  results 
in  the  equations: 

h(p)=r 

Dj  {  ElnBfx^n  i;Pi)~  h (p) }  =0  ,  l<j<m  where  h(p) 
is  the  reliability  function  of  the  coherent  system  and 
P*(P  »P  /  ,p ^  is  the  vector  of  component  reliabilities.  Myhre  and 
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Saunders  show  that  under  certain  conditions  p  is  a  decreasing  function 
of  <S  ( <5>0)  and  exists  as  a  fixed  point  of  a  contractive  map. 
Representing  this  fixed  point  as  p (5 ) ,  ln(L(p)  }=-x2(l)/2  may  be  turned 
into  an  equation  in  <5  .  As  ln{L(p(6))}  is  decreasing  in  6  ,  we  can 
bracket  -  X?l)/2  by  incrementing  from  zero.  Then  using  a  technique 
like  the  secant  method  to  solve  for  <5  we  may  find  the  confidence 
interval  [R,l]  ,  where  R=h(p(s))  . 

In  practice  ,  it  turns  out  that  the  contractive  map  may  be 
sensitive  to  changes  in  <S  near  the  solution  and  also  may  be  erratic 
when  the  map  is  started  far  from  the  fixed  point.  The  present 
subroutines  used  to  calculate  the  confidence  intervals  are  written  to 
extrapolate  (or  interpolate)  to  a  guess  at  the  new  fixed  point  each 
time  a  new  5  is  used.  By  doing  this  the  programs  have  run 
efficiently  and  trouble-free. 


The  Subroutines  for  General  Coherent  Systems 


There  are  two  subroutine  packages,  one  for  general  coherent 
systems  and  another  for  pseudo-series  systems  only.  Each  package  has 
three  subroutines,  LRSUB(  PLRSUB) ,  POSDEL,AND  CONMAP,  as  well  as  two 
functions,  XLAMBDA  and  RELFUN.  THEY  are  written  in  fortran  77  and  run 
on  the  CLAREMONT  COLLEGES  VAX  computer  when  combined  with  a  driving 
program  that  sets  up  the  data.  The  subroutines  in  the  package  for 
general  coherent  systems  and  RELFUN  share  a  common  block  of  data  that 
contains  the  structure  of  the  system,  the  test  data,  and  the  binomial 
component  test  data: 

Structure  of  the  System 
ORDER  is  the  number  of  distinct  components  in 
the  system 

ITYPE  is  a  parameter  of  the  type  of  system 
l:series 
2:parallel 
3:k  out  of  n 

Binomial  Test  Data 

WINS (I)  is  the  vector  of  number  of  successful 
trials  for  each  component 

TRIALS (I)  is  the  vector  of  number  of  times  each 
component  was  tested. 


Common  Computational  Variables 
DELTA  is  the  lagrange  multiplier 
OLDLAM  is  the  value  of  the  function  a( 
in  [M,S]  for  the  present  delta. 

I 


PDELTA(I)  is  the  vector  of  p^(<$)'s. 

PTILDE(I)  is  the  vector  of  maximum  likelihood 
estimates  of  the  p 

i 

LRSUB  sets  up  the  variables  and  calls  POSDEL  three  times  to 
compute  80, 85, and  90%  lower  confidence  bounds . POSDEL  solves  the 
equation  (p(  /2  by  incrementing  6  until  the  solution  is 

bracketed  and  then  interpolating  to  the  solution.  POSDEL  calls  CONMAP 
to  get  p(  6)  for  each  5  and  uses  XIAMBDA  to  compute  A(p(6))  once  p^S) 
is  known.  When  the  equation  is  solved,  control  returns  to  LRSUB  from 
POSDEL  and  tilSUB  computes  the  bound  by  plugging  the  solution  vector 
p(<5  )  into  RELFUN  ,  the  reliability  function  of  the  system. 


LRSUB 


LRSUB  is  the  subroutine  that  drives  the  programs  that  compute 
80, 85, and  90%  lower  confidence  bounds  for  general  coherent  systems 
using  the  likelihood  ratio  method.  LRSUB  sets  up  CHISQ(3)  as  x2(l) 
for  the  three  levels;  sets  up  p=.l  min{  WINS(I):  1<I<  ORDER  }  to  be 
used  as  an  increment  of  DELTA  in  POSDEL;  and  initializes  PDELTA(I)= 
PTILDE (I) =  WINS ( I ) /  TRIALS (I).  Then  LRSUB  calls  the  subroutine  POSDEL 
once  for  each  of  the  three  confidence  levels  with  the  appropriate 
q=-X2(1)/2;  computes  the  bound  by  using  RELFUN  and  the  returned 
PDELTA;  and  stores  the  computed  bounds  in  LCB(I) ,1=1,3.  LRSUB  calls 
POSDEL  ,  while  POSDEL  calls  CONMAP  and  XIAMBDA. 


POSDEL 


POSDEL  recieves  the  arguments  (P,Q)  and  uses  DELTA  and  the  other 
variables  in  the  common  block  to  solve  (p (6 ) ) =XLAMBDA (PDELTA) =Q= 
-X2(l)/2  as  an  equation  in  feDELTA,  where  PDELTA  is  a  vector  of  length 
ORDER  determined  by  CONMAP  for  each  DELTA  .  POSDEL  takes  into  account 
the  sensitivity  of  the  contractive  map  in  CONMAP  and  the  need  for 
computational  efficiancy  by  extrapolating  (or  interpolating)  from  the 
previous  values  of  PDELTA  to  get  a  new  PDELTA  which  corresponds  to  the 
new  DELTA  which  is  sent  to  CONMAP. 

To  keep  the  solution  process  clear  ,  POSDEL  is  written  in  two 
parts:  the  first  part  to  bracket  Q  by  finding  successive  DELTA  such 
that  XLAMBDA (DELTA-GAP) >Q>XLAMBDA (DELTA)  ;  and  the  second  part  uses 
the  secant  method  to  find  a  DELTA  such  that  ! XLAMBDA (DELTA)  -  QK.001. 
To  keep  track  of  the  past  results,  enabling  us  to  extrapolate  and 
interpolate  for  PDELTA  ,  LWBD(I)  ,  and  UPBD(I),  1=1, ORDER  are  set  up  to 
contain  the  PDELTAs  that  gave  the  closest  upper  and  lower  bounds  on  Q 
(LWBD(I)  is  initialized  as  0.0  and  UPBD(I)  as  x^/n^  for  the  first 
bound  and  as  the  solution  to  the  last  bound  for  the  second  and  third 
bounds) . 

In  part  one  DELTA  is  incremented  by  the  value  of  GAP  ,  which  is 
initially  the  value  of  P  ,sent  by  LRSUB,  and  for  each  new  DELTA  the 
vector  of  PDELTA' s  is  found  by  CONMAP  and  OLDLAM  is  computed  by 
XLAMBDA.  Part  one  iterates  through  values  of  DELTA  until  Q>OLDLAM 
( [M,S]  shows  that  XLAMBDA  is  a  decreasing  function  of  DELTA).  Once  we 
have  PDELTA  for  one  DELTA  ,  we  extrapolate  for  a  guess  at  the  next 
PDELTA  so  that  it  corresponds  to  the  new  DELTA  =DELTA  +GAP  ;  let 

p  =P  (5 1 ) ,  p  =P(  S 3=P(  <54-  £  and  assume  P(  3  is  linear  with  p.  =P(  $,■ fa  ) . 

12  4  *■ 

Then  Pj  -p2=slope*  &  =  (p  2~Pi )  /A  *A  and  ^  =2p  •  Thus  the  PDELTA  sent 


to  CONMAP  is  given  as 

PDELTA  { I )  =2*UPBD ( I )  -OLP  (I )  ,  KKORDER  where  OLP  is  the  PDELTA  from 
the  previous  step. 


If  CONMAP  fails  when  called  in  part  one  ,  GAP  is  halved  so  that 
POSDEL  will  approach  the  solution  region  more  slowly.  An  analysis 
similiar  to  the  above  shows  that  the  extrapolated  PDELTA  should  then 
be: 

PDELTA (I)  =3/2  UPBD(I)  -1/2  OLP  (I)  ,  KKORDER.  To  keep  LWBD  and  UPBD 
correct  they  are  updated  (after  saving  UPBD  in  OLP  )  before  completing 
the  loop.  Until  Q  is  bracketed  the  new  PDELTA  returned  by  CONMAP  is 
saved  in  UPBD  while  PDELTA  is  stored  in  LWBD  when  OLDLAM<Q  and  the 
program  moves  to  part  two. 

Upon  entering  part  two  ,  Q  is  bracketed  and  hence  we  may  estimate 
the  solution  by  the  interpolation: 

<5  =RR(52-5]+  where  RR=(Q-A1  )/(/  *  Using  the  same  interpolation 
P(  &)=  RR^-p^+p^  or  PDELTA (JJ)=  UPBD (JJ)+RR*  (LWBD (JJ) -UPBD (JJ)  , 
1<JJ<0RDER.  Once  DELTA  and  PDELTA  are  set  up  ,  CONMAP  is  called  and 
if  it  does  not  fail  LWBD  and  UPBD  are  updated  depending  on  whether  the 
new  OLDLAM  is  less  than  Q  or  larger  than  Q.  If  CONMAP  fails  in  part 
two  ,  RR  is  halved  and  CONMAP  tried  again  until  CONMAP  works  or  has 
failed  nine  times  (CONMAP  works  most  easily  if  it  does  not  have  to 
find  PDELTA  corresponding  to  a  large  change  in  DELTA) .  DELSAV  and  DL 
are  also  updated  before  a  new  interpolation  is  completed.  The  loop  in 
part  two  is  thus  implementing  the  secant  method  and  continues  until 
IQ-OLDLAMK.OOl  at  which  time  control  returns  to  LRSUB  where  the  bound 
is  computed. 
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CONMAP 


CONMAP  receives  as  arguments  FAIL, UPBD,  LWBD,  and  GAP,  and  uses 

DELTA  and  the  other  variables  in  the  common  block  to  compute  the  fixed 

point  of  the  contractive  map  A((ja  )  ,<5 )  of  Theorem  2-4  of  [M,S]  :  if 

(q  )=A((p  ),<$)  ,  and  if  a  (p)=p  *D  h(p) 

■**  i.  J  J  J 

then  q.  =  (x.  -a  (p)  )/(n  -a  (p) ) .  In  [B,P]  it  is  shown  that 

j  j  J  j  j.; 

D  h  (p)  =h  (p  q) -h  (p  x)  ,  where  Pj  is  the  vector  of  p  's  with  p  =1  and  pQ 
J  1  i  J 

is  the  same  vector  with  p  =0.  In  the  program  AJP=A  (p)  for  each  j  and 

j  J 

q  =A(J)= (WINS (J)-DELTA*AJP)/ (TRIALS (J) -DELTA* AJP) .  This  is  not 

j 

possible  if  WINS(J)  is  less  than  DELTA*AJP  and  in  that  case  A(J)  is 
left  at  the  previous  value  (which  is  PDELTA(J)  during  the  first  time 
through  CONMAP) .  If  the  new  value  of  A(J)  is  not  in  the  interval 
[UPBD ( J) , LWBD ( J) ]  it  is  set  equal  to  (PDELTA(J)+BD)/2  where  BD  is  the 
bound  it  exceeded. 

As  the  new  A(J)  are  computed  ,  DIST  is  updated  so  that  it  is  the 
distance  between  the  new  and  old  vector  values  of  the  contractive  map. 
The  map  fails  if  this  distance  doesn’t  contract  at  each  step  ,  which 
is  the  result  if  DIST>DIST2.  Otherwise  DIST2  is  updated  and  the 
process  continues  until  we  find  a  fixed  point  (DISTTOL=.0005)  and 
CONMAP  returns  successfully.  HOLD(I)  is  used  to  save  the  initial 
value  of  PDELTA  to  be  returned  if  the  contractive  map  fails.  If  GAP 
is  small  (GAP<. 00001) ,  PDELTA  is  determined  by  averaging  UPBD  and  LWBD 
since  little  change  can  be  expected  from  the  contractive  map  in  such  a 


RELFUN 


RELFUN  is  a  subfunction  that  computes  the  reliability  of  the 
given  type  of  system.  ITYPE  determines  the  type  of  system:  I TYPE  = 
If 2,  or  3  if  the  system  is  respectively  series,parallel  ,  or  k  out  of 
n  system.  This  function  takes  a  vector  of  probabilities,  PROB  ,  as 
its  argument  and  computes  the  reliability  using  the  appropriate 
formula: 

ITYPE=1:  series  system;  r=  nprob(i) 

ITYPE=2:  parallel  system:  r=  1-  II  (l-prob(i) ) 

ITYPE=3:  k  out  of  n  system  :  r=EE(i;n;p), 
where  E(i;n;p)  is  the  probability  of 
exactly  i  successes  out  of  n 


tv 

C******-**********  THIS  IS  THE  SUBROUTINE  LRSUB  ****■*•»*•»****■*■**•*** 
C  IT  RUNS  THE  LIKELIHOOD  RATIO  METHOD 

C  IT  DOES  NOT  ADJUST  THE  FAILURES 

C  IT  WORKS  FOR  ALL  SYSTEMS 

C  WHOSE  RELIABILITY  FUNCTION  CAN  BE  WRITTEN  IN  RELFUN 

C  NOW  ON  VAX.  2/23/30 

C  ALL  IS  WORKING  WELL  WITH  THE 

C  SUBROUTINES  POSDEL , CONMAP 

C*-**-#*  HH*************************-* **#**■**■*#*****■*#***■*********•*■*•#•■*•*■* 

SUBROUTINE  LRSUB (LCB) 

DIMENSION  CHI ISQ(3) 

INCLUDE  'lrcomr.for' 

INTEGER  ORDER 
REAL  LCB(3) 

DATA  CHI I SQ/. 709,  1.074,  1.642/ 

C****  P  GIVES  US  AN  INITIAL  INCREMENT  BY  WHICH  WE  CHANGE 
C***-*  DELTA  TO  LOCALIZE  ITS  VALUE 
P=1 . E15 

DO  40  1=1, ORDER 

IF ( WINS ( I  2 . LE. P)P=WINS( I ) 

PT I LDE ( I ) =W I NS  < I ) / TR I ALS ( I ) 

PBELTA< I ) =PTILD£( I ) 

40  CONTINUE 

P=P/10.0 

DELTA=P 
GLDLAM=0.0 
DO  50  1=1,  3 
Q=-1.0*CHII3Q< I)/2 
CALL  POSDEL  <  BOUND , P , Q ) 

BOUND=RELFUN  <  PDELTA ) 

LCB (IJ  =BOUND 
50  CONTINUE 

write <23, 1514) <  LCB ( I ) , 1  =  1,3) 

1514  format  (2x,  ''The  Bounds  are  s ', 3< f3.5> 2x) ) 

RETURN 
999  END 

C************************  END  OF  LRSUB  **■»****■»■»**■»■»***■»***•»**•» 

C 

C 

C 

C 

C  SUBROUTINE  POSDEL  ITERATES  THROUGH  VALUES  FOR  DELTA  PLUS 
C 

C*#**##***##*****#*#-**#*##  POSDEL  *****,‘********************* 

C  POSDEL  HAS  TWO  MAIN  parts!  in  part  1  the  value 

C  of  delta  is  incremented,  with  conmap  finding 

c  the  vector  of  p(i)s, until  p  has  been  bracketed, 

c  Once  a  vector  of  p<i)s  is  known  we  extrapolate 

C  to  set  a  new  vector  to  start  conmap  with(the  vector 

C  corresponds  to  the  new  value  of  delta).  This  makes 

C  conmap  less  likely  to  fail.  If  conmap  fails  the  size 

C  of  the  increment  for  delta  is  halved. 

C  Part  2  interpolates  to  get  a  new  delta  and  also 

C  interpolates  the  p(i)s  to  get  a  corresponding  vector 

C  to  start  conmap.  The  interpolation  is  made  possible 

C  by  storing  the  p<i)s  that  correspond  to  the  upper 

C  and  lower  bounds  on  Q. 

C  JBL  3/6/SO 

C  ******************************************************************* 
SUBROUTINE  POSDEL (BOUND, P, Q) 


»»!/'•  i  i»hp*  y  i  -  - 


p- / 


DIMENSION  UPBD (100),  OLP <100) 

REAL  LWBD(IOO) 

INTEGER  ORDER 
LOGICAL  FAIL 

C*****************-***  INITIALIZING  ***********■»■»**■»■»**■*****•»# 
GAP=P 

DO  121  1=1, ORDER 
LWBD< I )=0. 0 

UPBD ( I ) =PDELTA  < I ) 

121  CONTINUE 

FA I L=. FALSE. 

C***-*****#*******-****  PART  1  •*■»**■*■»*■»*****•»■»#■»**•»■»■»■»*■»*•»* 

C  START  OF  THE  SEARCH  FOR  A  SOLUTION  REGION  BY 
C  INCREMENTING  DELTA  BY  GAP  UNTIL  THE  SOLUTION  IS  BRACKETED 
C 

DO  20  J= 1,200 
OLAMB=OLDLAM 
DELSA V=DELT A 


DELTA=DELTA+GAP 

IF(J.GT.l)  THEN 

C  *•**•»-»•*  here  tue  extrapolate  for  pdelta  ******* 

DO  1042  JJ=1, ORDER 
PDELTA ( JJ ) =2*UPBD ( J  J ) -OLP ( JJ ) 

I F ( FA I L ) PDELTA ( J J ) = 1 . 5*UPBD ( J J ) - . 5*0LP ( JJ ) 
1042  CONTINUE 

END  IF 

CALL  CONMAP (FAIL, UPBD , LWBD , GAP ) 

1044  IF (FAIL)  THEN 

c  - block  1 - 

TYPE  9370 

9870  format ('  conmap  failed  ') 

DELT  A=DELT  A-GAP 
GAP=GAP/2. 0 

c  -  block  2 - 

ELSE 

OLDLAM® XLAMBD ( TR I ALS , PTILDE, PDELTA , ORDER ) 
d  rel =re 1 f un ( pdel ta ) 

d  write(5, 1  111 )  del ta, ol dl am, p, rel 

d  write(5, 1112)  (pdel ta( IT3T) , ITST=1 , order > 

1111  format (2x,4f 10.4, '##' , / , 3f6. 4) 

1112  f ormat ( f 8. 4 ) 


IF(ABS(OLDLAM-Q) .LT. .001 ) 


RETURN 


19 


DO  19  JJ=1, ORDER 
OLP ( JJ ) =UPBD ( J J ) 

IF ( OLDLAM. LT .  Q )  LWBD ( J J ) =PDELTA ( J J ) 
IF ( OLDLAM. GT . Q )  UPBD ( JJ ) =PDELTA ( JJ ) 
CONTINUE 

IF ( OLDLAM. GT. Q ) OLAMBU=OLBLAM 


c  - 

C  ****  THIS  IS  OUR  GOAL  FOR  THIS  PART 

IF(OLDLAM.LT.Q)  GOTO  30 

c - end  of  block  2 - 

END  IF 

20  CONTINUE 


WRITE(5, 2) 

2  format (4x, 'No  solution  after  100  steps  in  part  1 

BOUND®. 01 
RETURN 

C****************  END  OF  PART  1  *»■»****•»*****■»***•»■»#•»** 
C 


) 


C  THIS  IS  THE  START  OF  THE  SECTION  THAT  SOLVES  FOR 
C  DELTA  USING  THE  REGLILA  FALSI  METHOD 

C  AND  we  interpolate  for  the  Pdelta  (between  upbd 

C  and  1 wbd. ) . 

C********  *****************************  *********************** 


C 

30 


1313 


****  INITIALIZE  THE  VARIABLES  ***************** 
del tau=de1 ta-gap 
del tal =de 1 ta 
olamb 1 =o 1 dl am 

DO  40  J=l, 100 

****  here  we  interpolate  for  delta  and  pdelta  ***** 
IF ( .NOT. FAIL)  RR= ( Q-OLAMBU ) / ( OLAMBL-QLAMBU ) 

IF(FAIL)  RR=RR/2. 0 

DELT  A=DELTAU+RR* ( DELTAL-DELTAU ) 

DO  1313  JJ=1, ORDER 

PDELTA ( J J ) =UPBD <  J J ) +RR* <  LWBD ( J J ) -UPBD ( JJ ) ) 
CONTINUE 

GAPP=DELT  AL-DELT  AU 


9871 
c 


CALL  CONMAP (FAIL*  UPBD i LWBD*  GAPP ) 

IF ( .NOT. FAIL)  THEN 

-  block  1  - 

IT=0 

OLDLAM=XLAMBD<  TRIALS, PTILDE, PDELTA. ORDER ) 


IF ( (OLDLAM-Q) . 1 t. .001 ) RETURN 


33 

c 


39 

c 


IF(oldlam.lt.Q)  THEN 

—  block  a - 

DELT  AL=DELT  A 
OLAMBL=OLDLAM 
do  38  Jj=li order 
LWBD ( JJ ) =PDELTA <  J J ) 

CONTINUE 

ELSE 

-  block  b  - 

DELT  AU=DELT  A 
OLAMBU=OLDLAM 
DO  39  JJ=1 , ORDER 
UPBD ( J J ) =PDELTA ( J J ) 

CONTINUE 

END  IF 

****  CHECK  TO  SEE  IF  DELTAS  ARE  TOO  CLOSE  ***** 

IF( (del tal-del tau) . 1 t. IE-5)  then 
WR I TE (5*3) 

FORMAT  ( '  DELTA  CAN  NO  LONGER  BE  sisnif.  INCREASED.’') 

RETURN 

END  IF 


c 

40 


ELSE 


-CRISIS  ROUTE  - 

IT=IT+1 

IF( IT. st.9)stop 


END  IF 


CONTINUE 
WRITE(5, 1) 

FORMAT ( *  NO  SOLUTION  FOUND  AFTER  100  times  thru  part  2') 


BOUND®. 01 
RETURN 


c 


unsuccessful  exit 


c  THIS  COMPUTES  THE  CONTRACTIVE  MAP 

C  OF  EQN.  2.1  IN  MYHRE-SAUNDERS  PAPER. 

C  UPDATED  TO  JUST  KEEP  A<J>  BETWEEN  THE 

C  BOUNDS  ON  PDELTA.  JBL  2/16/80 

C  THIS  VERSION  IS  FOR  GENERAL  COHERENT 

C  SYSTEMS  AND  IS  USED  BY  LRSUB.  JBL  3/22/30 

C************************  CONMAP  *#**■***■*■**■*■*■#**■*■*•*•*•**•»*■*•***■*** 

SUBROUT I NE  CONMAP < FAIL, UPBD, LWBD » GAP ) 

INCLUDE  "  1  rcomr .  FOR’' 

INTEGER  ORDER 
LOGICAL  FAIL 

DIMENSION  HOLD ( 1 00 ) , A  < 1 00 ) , UPBD < 100) 

REAL  LWBD! 100) 

FAIL®. FALSE. 

TOL=. 0005 

C  ****  IF  GAP  IS  VERY  SMALL  WE  SIMPLY  AVERAGE  THE  BOUNDS  OF  PDELTA. 

IF ( GAP. LT. IE-4 ) THEN 
DO  999  J=l, ORDER 
PDELTA ( J )  =  ( LWBD <  J ) +UPBD <  J ) ) /2 . 0 
999  CONTINUE 

RETURN 

c  - 

ELSE 

C  MAIN  PART  OF  WORK. 

DO  10  J=l, ORDER 
HOLD  <  J ) =PDELT A  <  J ) 

A  <  J ) =HOLD ( J ) 

10  CONTINUE 

C*****  SOLVE  THE  CONTRACTIVE  MAP 
DIST2=1 . 0 

100  DO  40  I CE= 1,1000 

DIST=0. 

c  -  compute  delaj  for  each  component  - 

DO  30  J=l, ORDER 

IF<WINS< J) .EQ. TRIALS! J) >A( J)=l. 

IF  <A< J) .NE. 1.0)  THEN 
temp=pdel ta (J ) 
pde 1 ta  < J )=0. 0 
p0=re1 f un  <  pde 1 ta ) 
pdel ta< j )=1.0 
pl=re1fun(pdelta) 
pde 1 ta  <  j )  =  temp 

de 1 aJ=de! ta*pde 1 ta  <  j ) *  <  pI-pO ) 

C  ***  IF  DELAJ>=WIN3,  DO  NOT  ASSIGN  A(J>,  IT  WILL  BE  NEG. 

IF< DELAJ.LT. WINS! J) )  THEN 

A! J)=! WINS! J) -DELAJ) /(TRIALS ! J) -DELAJ) 

C  ***  WE  MUST  MAKE  SURE  A(J)  STAYS  BETWEEN  THE  BOUND 

IF!A! J) .LT.LWBD! J) )A( J)=  (PDELTA! J)+LWBD( J) ) /2. 0 
IF! A ! J) . GT. UPBD ! J ) ) A < J )=! PDELTA ! J) +UPBD! J ) ) /2. 0 
DIST=DIST+ABS(A( J) -PDELTA! J) ) 

c  - 

END  IF 

END  IF 


30  CONTINUE 

c  - - - - new  pdel  t a s  are  done 

DO  70  L«l. ORDER 
PDELT  A  <  L ) ® A ! L ) 

70  CONTINUE 


*****  IF  APPROXIMATIONS  HAVE  CONVERGED  ENOUGH,  WE  STOP  ITERATING 
DIST*D I ST/ORDER 


o  « 


c 

c 


41 

c 


50 


40 


IF(DIST.LT.TOL)  RETURN 

*****  IF  DI3T  IS  INCREASING,  CONMAP  FAILS  AND  WE  GO  TRY 
*****  ANOTHER  DELTA  AND  KEEP  THE  OLD  PDELTA'S. 

IF (DIST. LT. DIST2 )  THEN 
DIST2=DIST 


ELSE 

DO  50  J=1 ,  ORDER 
PDELT  A  ( J ) =HOLD ( J ) 
FAIL=. TRUE. 

RETURN 


END  IF 

CONTINUE 
END  IF 
END 


*********************************************************** 


Subroutines  for  Pseudo-series  Systems 


The  subroutine  package  for  pseudo-series  systems  consists  of 
three  subroutines,  PLRSUB,POSDEL,  and  CONMAP  as  well  as  two  functions, 
XLAMBDA  and  RELFUN.  These  subroutines  perform  the  same  tasks  as  those 
in  the  other  package  ,  but  are  specialized  for  pseudo-series  since  the 
reliability  function  of  a  psuedo-series  is  given  as  h(p)=£  \jllp  J 
where  a ^>0  if  the  i-th  component  is  in  the  j-th  subpopulation  and 
a^.=0  if  not.  See  [IDS  1980]  for  a  derivation  of  this  function  that 
shows  how  it  models  some  missile  systems.  From  this  form  of  the 
reliability  function  it  follows  that  D^i(p)*p^=  £  w.altpa.  In  this 
package  the  part  of  the  common  block  describing  the  structure  of  the 
system  is  expanded  so  that  it  contains: 

ORDER  is  the  number  of  distinct  components  in 


the  system. 

LLL  is  the  number  of  subpopulations  in  the  system. 

W (J)  is  the  vector  of  subpopulation  weights. 

AL(I)  is  the  vector  of  component  exponents. 

IG (I, J)  is  the  matrix  of  indices  for  the  inclusion 
of  component  j  in  the  i-th  subpopulation.  The  rest  of  the 
common  block  is  the  same  as  in  the  other  package. 


PLRSUB  performs  exactly  the  same  role  in  this  package  as  LRSUB 

does  in  the  other  package.  POSDEL  is  the  same  as  the  other  POSDEL, 

while  CONMAP  is  almost  the  same  except  for  the  lines  computing  AJP. 

Since  a  (p)=p  *D  h(p)=  £  wodipa  and  is  not  p  (h(a  )-h(p  )),  AJP  is 
k  k  k  j  0  1 

computed  directly  in  this  version  of  CONMAP.  XLMBDA  is  the  same  as 
before,  but  of  course  RELFUN  is  the  implimentation  of  the 
pseudo-series  reliability  with  no  possibility  of  modelling  anything 
else. 


m 


tv  pj  ust  1 r. for 

C**************  THIS  IS  PJUSTLR  **************** 

C 

C  REVISED  2/8/80  JBL 

C*************************************************** 

REAL  LCB(3) 

INTEGER  ORDER 
INCLUDE  ' 1 rcomp. FOR' 

CALL  PGETDAT (WINS) 

CALL  PLRSUB  <  LCB ) 

999  END 

C  ************************************************* 

C 

c 

c 

C******************  SUBROUTINE  PGETDAT  ***************** 

C  ******************* 

C  THIS  WILL  SET  UP  THE  DATA  FOR  COMPUTING  LIKELIHOOD 

C  RATIO  LOWER  CONFIDENCE  BOUNDS  FOR  PSEUDO-SERIES  BY 

C  ASKING  FOR  THE  INPUT  AND  OUTPUT  FILES  AND  OPENING  THE  FILES. 

C  THE  INPUT  FILE  MUST  CONTAIN  IN  ORDER: 

C  NUMBER  OF  UNIQUE  COMPONENTS  (N)  AND  SUBSYSTEMS  <L) 

C  WEIGHTS  FOR  THE  L  SUBSYSTEMS 

C  EXPONENTS  FOR  EACH  UNIQUE  COMPONENT 

C  NXL  MATRIX  INDICATING  WHICH  OF  THE 

C  L  SUBSYSTEMS  EACH  COMPONENT  IS 

C  COMMON  TO. 

C  SUCCESSES  AND  TRIALS  FOR  COMPONENTS 

C  1,2,. ..,N 

C 

C******************************************************************** 
SUBROUTINE  PGETDAT <0W I NS) 

DIMENSION  OWINS(IOO) 

INTEGER  ORDER 
DOUBLE  PRECISION  FILNAM 
INCLUDE  '1 rcomp. FOR' 

15  FORMAT (2i2) 

TYPE  1313 

1313  FORMAT ( '  WHAT  FILE  HAS  THE  SYS.  DATA  ',*) 
accept  1314,  FILNAM 

1314  FORMAT ( A10) 

OPEN < UNI T=22, name=f i 1 nam, ACCESS-' sequential ' , tvpe='ol d ' ) 

TYPE  1315 

1315  FORMAT <  '  WHAT  FILE  SHOULD  THE  OUTPUT 
2  BE  FILED  IN?'*) 

accept  1314, FILNAM 

open  C  unit=23,  name=f i  1  nam,  access3'  sequential  '  ,  tYPe='new'  ) 

READ  <22,15) ORDER , LLL 
IF  (ORDER  .EQ.  0)  GOTO  30 
READ (22,413) <W( I) , 1=1, LLL) 

413  FORMAT ( 12F) 

write (23, 1511 ) 

1511  format (2x, 'THE  FOLLOWING  IS  RUN  ON  THE  VAX,',/ 

1  , 2x» 'usins  JBEPOS, JBCONM  .') 

call  IDATE( IDl»id2,id3> 
write(23, 1512)  idl,id2,id3 

1512  format(2x, 'The  date  is  ', i3, 'i '» i3, '«', i3> 
write (23, 1513)0RDER, LLL 

1513  format (2x, /, 15x, 'THE  SYSTEM  DATA  IS  :',/, 

1  4x» 'THE  order  is  '»i3, '  the  number  of  sub  is  ',13) 

WRITE ( 23, 512) 

512  FORMAT < / , '  THE  WEIGHTS  ARE:  ') 

WRITE(23, 413) ( W< J) , J=1 , LLL ) 

_ _  PEAD<22, 414)/ AL.<  J>,  J=1 ,  ORDER). _ _ _ 


414 


FORMAT ( 1 OF) 

WRITE ( 23,  515 ) 

515  FORMAT (/,'  THE  INDICES  THE  EXPONENTS 

DO  119  J=l, ORDER 
READ (22, 415) ( IG< I , J) , I=1,LLL) 

WRITE (23, 416)J, ( IG ( I , J) , 1=1 , 12) , AL( J) 

415  FORMAT ( 131 ) 

416  FORMAT (1313, FIO. 6) 

119  CONTINUE 

DO  30  1=1, ORDER 

READ(22, 125)  OWINS( I ) , TRI ALS ( I ) 

30  CONTINUE 

125  FORMAT (2F) 

WRITE (23, 1225)  ( OWINS ( I ) , 1=1, ORDER) 

1225  FORMAT ( *  WINS:  '.15F7.3) 

WRITE ( 23, 1226)  (TRIALS(I) . 1=1, ORDER) 

1226  FORMAT ( "  TRIALS:  y15F7. 1) 

RETURN 

END 

C******-*****************  XLAMBDA  ******■»********•*■**•*■#•**■**•»■**•***■#••»■ 

C  THIS  SOLVES  EQUATION  1.4 

C  TO  FIND  LAMBDA (DELTA) 

C 

FUNCTION  XLAMBD( TRIALS, PTILDE, PDELTA, I ) 

DIMENSION  TR I ALS ( 1 00 ) , PT I LDE ( 1 00 ) , PDELT A ( 1 00 ) 

S=0. 

DO  10  J=1 , I 

•  IF(PTILDE(J).EQ. 1.0)GOTO  10 

2L0G1=AL0G(PDELTA( J) /PTILDE ( J) ) 

ZL0G2=AL0G ( ( 1 . -PDELTA ( J ) ) / ( 1 . -PTILDE ( J) ) ) 

S=S+TRI ALS ( J) * ( PTILDE ( J ) *ZLOG 1 + ( 1 . -PT I LDE ( J ) ) *ZL0G2 ) 

10  CONTINUE 

XLAMBD=S 
RETURN 
END 

C*****************-************* *******************************  ■**■*■** 

C 

c 

c 

c 

C #*###*#####*##*#####**  RELFUN  ********•*■***•****■*•*■*-*#**■*•■*■*■»*■#••*••**** 

C  THIS  IS  THE  REL.  FUNCTION  FOR  PSEUDO-SERIES 

C 

C 

FUNCTION  RELFUN (PROB) 

INCLUDE  ' 1 rcomp. FOR' 

DIMENSION  PROB ( 1 00 ) 

INTEGER  ORDER 
S1=0.0 

DO  100  L=1 , LLL 
Pl-1.0 

DO  90  J-l, ORDER 

P1=P1*PR0B( J)**(AL< J)*IG(L,  J) ) 

90  CONTINUE 

S1=S1+W<L)*P1 
100  CONTINUE 

RELFUN=S1 
RETURN 

“ .  END 

'  *  * 


tv  Plrsub.for 

C****************  THIS  IS  THE  SUBROUTINE  PLRSUB  ****************** 
C  IT  RUNS  THE  LIKELIHOOD  RATIO  METHOD 

C  IT  DOES  NOT  ADJUST  THE  FAILURES 

C  IT  WORKS  ONLY  FOR  PSEUDO-SERIES. 

C  ALL  IS  WORKING  WELL  WITH  THE 

C  SUBROUTINES  POSDEL, CONMAP 

C  WORKING  ON  THE  VAX,  5/9/SO 

C****************************************************************** 

SUBROUTINE  PLRSUB <LCB) 

DIMENSION  CH 1 1 SQ ( 3 > 

INCLUDE  ' 1 rcomp. for' 

INTEGER  ORDER 
REAL  LCB ( 3 ) 

DATA  CHI ISO/. 709,  1.074,  1.642/ 

C****  P  GIVES  US  AN  INITIAL  INCREMENT  BY  WHICH  WE  CHANGE 
C****  DELTA  TO  LOCALIZE  ITS  VALUE 
P=1 . E15 

DO  40  1=1, ORDER 

IF (WINS ( I ) . LE. P )P=WINS ( I ) 

PTILDE ( I ) =WINS ( I ) /TRIALS  < I ) 

PDELTA( I )=PTILD£( I) 

40  CONTINUE 

P=P/10.0 

DELTA=P 
0LDLAM=0. 0 
DO  50  1=1,  3 
Q=— 1 . 0*CH IISQ(I)/2 
CALL  POSDEL ( BOUND, P,Q) 

BOUND=RELFUN ( PDELTA ) 

LCB ( I ) =BOUND 
50  CONTINUE 

RETURN 
999  END 

C************************  END  OF  PLRSUB  ************************ 

C 

C 

C 

C 

C*************************  POSDEL  *************************** 

C  POSDEL  HAS  TWO  MAIN  parts:  in  part  1  the  value 

C  of  delta  is  incremented,  with  conmap  finding 

c  the  vector  of  p(i)s, until  s  has  been  bracketed, 

c  Once  a  vector  of  p(i)s  is  known  we  extrapolate 

C  to  get  a  new  vector  to  start  conmap  with(the  vector 

C  corresponds  to  the  new  value  of  delta!.  This  makes 

C  conmap  less  likely  to  fail.  If  conmap  fails  the  size 

C  of  the  increment  for  delta  is  halved. 

C  Part  2  interpolates  to  9et  a  new  delta  and  also 

C  interpolates  the  p(i)s  to  set  a  corresponding  vector 

C  to  start  conmap.  The  interpolation  is  made  possible 

C  by  storing  the  p(i)s  that  correspond  to  the  upper 

C  and  lower  bounds  on  Q. 

C  JBL  3/6/SO 

C  NOW  WORKING  ON  VAX  2/28/SO 

C  ************************************************************ 

SUBROUT I NE  POSDEL ( BOUND ,P,Q> 

INCLUDE  ' 1 rcomp. FOR' 

DIMENSION  UPBD<100),0LP(100) 

REAL  LWBD(IOO) 

INTEGER  ORDER 

LOGICAL  FAIL  _ _ _ 


C********************  INITIALIZING  ************************** 
GAP=P 

DO  121  1=1, ORDER 
LWBD< I )=0. 0 

UPBD( I ) =PDELTA  < I ) 

121  CONTINUE 

FAIL=. FALSE. 

C**##*#»##*##*******#  PART  1  **************************** 

C  START  OF  THE  SEARCH  FOR  A  SOLUTION  REGION  BY 
C  INCREMENTING  DELTA  BY  GAP  UNTIL  THE  SOLUTION  IS  BRACKETED 
C 

DO  20  J= 1,200 
OLAMB=OLDLAM 
DELSAV=DELT  A 
DELT  A=DELTA+GAP 

IF  ( J. GT. 1 )  THEN 

C  ******  here  we  extrapolate  for  Pdelta  ******* 

DO  1042  JJ=1, ORDER 
PDELTA  <  J J ) =2*UPBD ( JJ ) -OLP ( J J ) 

IF ( FA I L ) PDELTA ( JJ) =1 . 5*UPBD  <  J J ) - . 5*QLP  <  J J ) 

1042  CONTINUE 

END  IF 

CALL  CONMAP (FAIL, UPBD , LWBD , GAP ) 

1044  IF (FAIL)  THEN 

c  - block  1 - 

TYPE  9370 

9370  formats  conmap  failed  ') 

DELT  A=DELT  A-GAP 
GAP=GAP/2.0 

C  -  block  2 - 

ELSE 

OLDLAM=XL AMBD (TRIALS, PT I LDE , PDELTA, ORDER) 

I F ( AB3 ( OLDLAM-Q ) . LT. .001 )  RETURN 

c  ++++++++++++++++++++++++++++++++++++++++++++++++-f-++++- 

DO  19  JJ=1, ORDER 
OLP ( J J ) =UPBD ( J J ) 

IF ( OLDLAM. LT.  Q  >  LWBD ( JJ ) =PDELTA ( JJ ) 

I F ( OLDLAM . GT . Q )  UPBD ( J J ) =PDELTA ( J J ) 

19  CONTINUE 

IF ( OLDLAM . GT . Q ) OLAMBU=OLDLAM 

c  - 

C  ****  THIS  IS  OUR  GOAL  FOR  THIS  PART 

IF (OLDLAM. LT. Q)  GOTO  30 

c - end  of  block  2 - 

END  IF 

20  CONTINUE 

c - 

WRITE<5,2) 

2  "  format (4x,  'No  solution  after  100  steps  in  part  1.*') 

BOUND=. 01 
RETURN 

C****************  END  OF  PART  1  *********************** 

C 

C******************  PART2  ********************************** 

C  THIS  IS  THE  START  OF  THE  SECTION  THAT  SOLVES  FOR 
C  DELTA  USING  THE  REGULA  FALSI  METHOD 

C  AND  we  interpolate  for  the  pdelta  (between  upbd 

C  and  1 wbd. ) . 

C************************************************************ 

C  ****  INITIALIZE  THE  VARIABLES  ***************** 

30  del tau=del ta-eap 

del ta1=del ta 
o1amb1=o1 dlam 


OOOOn  n  n  ft  ftft 


c  ****  here  we  interpolate  for  delta  and  pdelta  ***** 

I F  < . NOT . FA I L )  RR= ( Q-OLAMBU ) / ( OLAMBL-OL AMBU ) 

IF (FAIL)  RR=RR/2. 0 
DELTA=D£LTAU+RR*  ( DELTAL-DELTAU ) 

DO  1313  JJ=1, ORDER 

P DELTA <  JJ ) =UPBD <  J J ) +RR* ( LWBD ( JJ ) -UPBD ( JJ )  ) 
1313  CONTINUE 

GAPP=DELT  AL-DELT  AU 


CALL  CONMAP (FAIL, UPBD , LWBD , GAPP ) 


IF(. NOT. FAIL)  THEN 

c  - block  1  - 

IT=0 

9871  OLDLAM=XLAMBD( TRIALS, PT I LDE, PDELTA, ORDER) 


IF( (OLDLAM-Q) . 1 t. . 001 )RETURN 


c 

c 


IF(ol dl am. 1 1.  Q)  THEN 

—  block  a - 

DELTAL=D£LTA 
OLAMBL=OLDLAM 
do  38  Jj=l, order 


LWBD ( JJ ) =PDELTA ( JJ ) 
33  CONTINUE 


ELSE 

c  - block  b - 

DELTAU=DELTA 
GLAMBU=OLBLAM 
DO  39  JJ=1, ORDER 
UPBD ( JJ ) =PDELTA ( J J ) 

39  CONTINUE 

c  - 

END  IF 

C  ****  CHECK  TO  SEE  IF  DELTAS  ARE  TOO  CLOSE  ***** 

IF( (del tal-del tau) . It. IE-5)  then 
type  3 

3  FORMAT ( '  DELTA  CAN  NO  LONGER  BE  sisnif.  INCREASED. " 

RETURN 
END  IF 


ELSE 

c  - CRISIS  ROUTE - 

IT=IT+1 

IF( IT. 5t. 9) stop 

END  IF 

c  "  - 

40  CONTINUE 

WRITE (5, 1 ) 

1  FORMAT ( '  NO  SOLUTION  FOUND  AFTER  100  times  thru  part  2') 

BOUND=».  01 
RETURN 

c  -  unsuccessful  exit  - 

END 


C***********************  END  OF  POSDEL  *#»*«#«««******#*«*#**** 


) 


■♦••►-M-+++++4'++-t-+++++++-*-++-f  CONMAP 


►+++++++++++++ 


THIS  COMPUTES  THE  CONTRACTIVE  MAP 
OF  EQN.  2.1  IN  MYHRE-SAUNDERS  PAPER. 

THIS  IS  ONLY  GOOD  FOR  PSEUDO-SERIES  SYSTEMS 

JBL  10/31/73 

_ _ UPDATED „T0  JUST  ,KEEP_A(J)  BETWEEN  THE  . . 


C  BOUNS  ON  PDELTA.  JBL  2/16/80 

C************************  CONMAP  ****************************** 


C 


999 

c 

C 


10 

C ***** 
100 
c 


26 

25 

C 


C 


c 


30 

c 


70 

c 

.  c 

c 

c 

c 


i 


SUBROUT I NE  CONMAP (FAIL, UPBD , LWBD , GAP ) 

INCLUDE  rcomp.FOR-' 

INTEGER  ORDER 
LOGICAL  FAIL 

DIMENSION  HOLD< 100) , A ( 100 ) , UPBD ( 100 ) 

REAL  LWBD (100) 

FA I L=. FALSE. 

TOL=. 0005 

****  IF  GAP  IS  VERY  SMALL  WE  SIMPLY  AVERAGE  THE  BOUNDS  OF  PDELTA. 

IF(GAP.LT. IE-4 ) THEN 

DO  999  J=l, ORDER 

PDELT A ( J )  =  <  LWBD ( J ) +UPBD ( J ) ) / 2 .  0 

CONTINUE 

RETURN 


ELSE 

MAIN  PART  OF  WORK. 

DO  10  J=l, ORDER 
HOLD ( J ) =PDELT A ( J ) 

A( J)=HOLD( J) 

CONTINUE 

SOLVE  THE  CONTRACTIVE  MAP 
DIST2=1 . 0 
DO  40  ICE=1 , 1000 
DIST=0. 

- compute  delaJ  for  each  component - 

DO  30  J=l, ORDER 

IF ( WINS ( J) .EQ. TRIALS (J) >A(J)=1. 

IF  (A(J).NE. 1.0)  THEN 
AJP=0. 0 
DO  25  L=1 , LLL 
P 1=1.0 

DO  26  JL=1, ORDER 

P1=P 1 *PDELTA <  JL ) ** ( AL ( JL ) * I G ( L , JL ) ) 

CONTINUE 

AJP=AJP+W(L)*AL< J)*IG(L, J)*P1 
CONTINUE 
DEL A J=DELT A* A JP 

***  IF  DELAU>=W I NS ,  DO  NOT  ASSIGN  A(J),  IT  WILL  BE  NEG. 

IF ( BELAJ. LT. WINS ( J) )  THEN 
A( J) = ( WINS ( J)— DELAJ) / ( TRIALS ( J) -DELAJ) 

***  WE  MUST  MAKE  SURE  A(J) .STAYS  BETWEEN  THE  BOUNDS 
IF ( A ( J ) . LT. LWBD ( J ) )A(J)=  ( PDELTA (J)+LWBD(J) )/2.0 

IF(A(J). GT. UPBD( J) ) A ( J ) = ( PDELT A ( J ) +UPBD ( J ) )/2.0 
D I ST=D I ST+ABS  <  A ( J ) -PDELT A ( J ) ) 


END  IF 

END  IF 
CONTINUE 

-  new  pdeltas  are  done  - 

DO  70  L=l, ORDER 
PDELTA(L)=A(L) 

CONTINUE 

*****  IF  APPROXIMATIONS  HAVE  CONVERGED  ENOUGH,  WE  STOP  ITERATING. 
D I ST*D I ST /ORDER 

+++++++++++++++++++++++-M-+++++++++++++++++++++++++++++ 

IF(DIST.LT.TOL)  RETURN 

*****  IF  DIST  IS  INCREASING,  CONMAP  FAILS  AND  WE  GO  TRY 
*****  ANOTHER  DELTA  AND  KEEP  THE  OLD  PDELTA 'S. 

IF<DIST.LT.DI ST2 )  THEN 
DIST2-DIST 


41 


X 


ELSE 

DO  50  J=l, ORDER 

50  PDELTA  <  J ) =HCLD ( J ) 

FAIL=. TRUE. 

RETURN 


END  IF 

40  CONTINUE 

END  IF 
END 

C*-****-* **********■#•■*■#*-#•**-**■*■#••#•*■****■*■•»■*  ■****■»*******■»*#**■«•*■»*■»* 
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